{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c3a415eb20>]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAdpElEQVR4nO3deXxVZ73v8c8vOzOEOYFCBqYABToAgbZYO9lW2tPxaI/UDtapVk+vVq9DfTnV43BOPbde73lZ5aK3Wl/aYtVasaKdFU/HhBYoIUDDmBAIgcxz9t7P/SMbTEIou5Bk7bX29/165cVez1rZ+/ck5cvTZw2POecQERH/S/G6ABERGRoKdBGRgFCgi4gEhAJdRCQgFOgiIgGR6tUHT5o0yU2fPt2rjxcR8aUNGzYcds7lDrbPs0CfPn06ZWVlXn28iIgvmdneE+3TlIuISEAo0EVEAkKBLiISEAp0EZGAUKCLiASEAl1EJCDiCnQzW2Fm282s0szuHWT/JWbWZGYbY19fH/pSRUTk7Zz0OnQzCwEPAlcA1UCpma11zm0dcOjfnXPXDEONIiK+1NDWzaOl++jsjvRrL5k+gYvmDHpv0GmJ58aiZUClc24XgJmtAa4HBga6iEjSi0QdG6saqapv5zvrKqhr6cKs/zF3XTzLs0CfBlT12a4GzhvkuAvMbBNQA3zeOVc+8AAzuxO4E6CwsPCdVysiMsLCkSi7DrcRdQ7n4M39Tew53HZs//7GDsr2NHB0saD2ngiN7T0AzJuSw8MfXsb8qWNGpNZ4At0GaRu4zNHrQJFzrtXMrgaeAIqP+ybnVgOrAUpKSrRUkogkrM6eCL98ZS8P/fduapo6++1LMQil9EZjakoKF8/JJSfzH3G6pGg8s/JGc3b+WDJSQyNWczyBXg0U9NnOp3cUfoxzrrnP63Vm9iMzm+ScOzw0ZYqIDK3qhnZ21bXx/LZD1LV29du3v6GDjVWNABTnjebfrl9A7ugMAMZlp3P+zAnYwHmUBBBPoJcCxWY2A9gPrAQ+2PcAM5sC1DrnnJkto/fqmSNDXayIyKlo6exh75F2Xtp5mJ6I4/lth9hY1Ugk6gilGEUTs/tNRaSFUrjr4llMn5jN+5bkkxbyxxXeJw1051zYzO4GngJCwEPOuXIzuyu2fxXwfuCTZhYGOoCVTqtPi8gI236whV11rce2a5o6ea6ilpd29h9fThuXxW3nF3H5mZOZkTuKaeOyRrrUYWFe5W5JSYnT43NF5FS1doX5xct7aIqdgKxr7eL3b+xnYKRNHZvJlQumMDtvNJfOy2N8dhqZqSFSUhJvyiQeZrbBOVcy2D7PnocuIvJ26lq6eK6ilhe2HyISPX7gWXGghf2NHWSm9U6HhMy49bwibl5WSEpshiQtlMLMSaMScr57OCjQRWTE1bd1E45Gj2vviTieKT/IgaZOHnl1Hy1dYfJyMsjNyTju2PzxWfxg5bksnT5hJEr2BQW6iJy2SNTREzk+oAGizvHa7nrW7zjM3iNtHG7tYlN100nf88LZk7j7stksKRrvm5OSXlOgi0hcDrV0UtfSxZrXqqhuaD/WHnXwxr4GmjvDb/v9oRSjOG806akpfOqSWUw9wYnIM8/IYUmRRt2nQoEuIsc50NTBa7vrcQ4ONnfy7NZayvY2AJCemsLcyTn9bmdfPmsSZxeMxQa9DxEKJ2Rz3swJTBp9/NSJDB0Fuogc89fth/jfz77FpthNNUcVTczm4++ewbwpY3j3nEnk5WR6VKG8HQW6SJJrbO/mq09sYXN1E/vq2ymckM2dF83kqoVTGJedTmqKkT8+K2muFPEzBbpIgIUjUWoaO3nrUAu7+zxQau+Rdta/VUfUOZrae+jsiXLZvDzetzifT14yi/RUnYT0IwW6SAB0h6P8dkM1b+5vpK7lH88l2V7bQlV9x6Df8+7iSeSOziCUYtxyfhHnFowbqXJlmCjQRXzMOcfq9bu4/y/biDrITg8xfeKoYycs83IyuWP5DArGZ7FsxoRjTwgMpRjZ6frrHzT6jYr4RDgSpTsSpeJAM3/bXkdLV5jfbqimpTPMu4snccfy6Vw6N8+3t7TL6VOgiySozp4Ir+w6wnMVh+jsifD01lqaOnqO7U8xOKdgHDcvLeSmknydtBQFuohXusIRHiurpq75H4sn7G/sZO+R3pOXuw+3caStm/TUFMZkpjFvSg6XzM0jKy2FGxfnMzYrzavSJUEp0EWGSTTqONLW3a8tEnWs31FHfXs3z1XUUrqn92ado4Pr7LQQZ+WPJZRiLC4az3sXTOG9CyaTk6nwlpNToIsMocOtXfxuQzUv7TzCxqrGflMkA6WnpvC/bjqH9y/JH8EKJcgU6CJxeGZrbb+FE07k92/sZ9vBFiaMSufyMyezYOoY0gZc0z1r0igWFY4nlGK63luGlAJdJCYSdTxWVkVtc/8Fgfc3dPCbDdVxvUdGagqrbl3MioVnDEeJIm9LgS5J7fV9Ddz7u83Ut/UQjkZpbD9+iiQ1xVi5tICv/NOZx67jPpFQio3oKu8ifSnQJXA6eyLsPMH0SMWBFn74/Ft0h3uf3V3X2sWUsZlcuWAyAOcWjOOmQea0dUmg+IECXQKhtStMVX07L1Ye5qd/383BAdMmfc2dnHNslZsxWWn866WzmTAqfaRKFRk2CnTxtdrmTr7zpwpe2H6IltgCC9PGZfG995896HXaITOWz56o294lkPRftfhS5aEWvv/MDjbua6S+vZuLinNZsXAKM2JXkIgkIwW6+NLDL+3lma21nDVtLD+8ZTGLFeIiCnTxn/buMM9W1HLJ3Dx+cnuJ1+WIJAzd1SC+4ZzjsdIqrvvhi9Q2d3LzsgKvSxJJKBqhiy845/j3P29j9fpdTByVzqpbl3DZvMlelyWSUBTo4gvff2YHq9fv4vYLirjv2gV65rfIIBToktB6IlH+86ntrF6/ixsXTeOb1y3QTT4iJ6BAl4TU0tnDn7cc5OnyWp6tqGXhtDF8acU8hbnI21CgS8Kob+vm87/ZxKaqRlo6w3RHem/P/9KKeXzyklkeVyeS+BTokhD2HG7js49tpLymmWvPnkp2eoj3nJnHuQXjGJet2/JF4qFAF88daunk5p+8wpG2bv7jn8/inxdrwQeRU6FAF0855/j0o2/Q2N7D459czsJpY70uScS3FOjimR/9tZLfbahmZ10b375hocJc5DTpTlEZcc45vruugu/9ZTvjs9P5xMUzWblUd32KnC6N0GVEbaxq5J41b7DnSDu3nV/EfdctOOkqQCISHwW6jJja5k4+8vNSMlJT+OKKudx10Szd8SkyhOIKdDNbAfwfIAT81Dn3Hyc4binwCvAB59xvh6xK8aXucJSqhnZaO8M8+EIlWw8009Ed4bFPXMDsvNFelycSOCcNdDMLAQ8CVwDVQKmZrXXObR3kuPuBp4ajUPGP0j31fPOP5exv6KAhtujy6IxUzs4fy33XLlCYiwyTeEboy4BK59wuADNbA1wPbB1w3P8AfgcsHdIKxTdqGju4b205G/Y2kBoyLpqTy/JZE8lMC7GoYDyFE7O9LlEk0OIJ9GlAVZ/tauC8vgeY2TTgRuAyFOhJ5+cv7ubXZdUcae2ipTPMkqLxfO7KOVpFSGSExRPog521cgO2fwB8yTkXebuHJ5nZncCdAIWFhfHWKAnqT5sP8MDT29lzpI0zzxjDOQXjuHlZgZ5TLuKReAK9Guh7kXA+UDPgmBJgTSzMJwFXm1nYOfdE34Occ6uB1QAlJSUD/1EQH6lr6eLexzczeUwmdyyfwWevKCYnM83rskSSWjyBXgoUm9kMYD+wEvhg3wOcczOOvjaznwNPDgxzCQbnHI++VsWvS/fR1RNl9W1LmJmrk5wiieCkge6cC5vZ3fRevRICHnLOlZvZXbH9q4a5RkkAFQea2VTVyAPP7KCupYsxmal858aFCnORBBLXdejOuXXAugFtgwa5c+6O0y9LEoVzjrWbarjn1xtxDqaOzeTbNyzklvMKtdiESILRnaJyQjvrWvn2k1t5YXsdE0el8/MPL6N48mgy00JelyYig1Cgy6B+/Ned3P+XbQB8oKSAL6yYy6TRGR5XJSJvR4Eux/xk/S6efPMAAFv2NzE7bzSrbl3M7LwcjysTkXgo0JPcC9sP8cir+/jbjjq6w1EWThvDxFEZXH/OVL56zXwmjNLybyJ+oUBPQpurG/nh85U8W1FL1EEoxbhx0TTmTcnhtguKyEjVHLmIHynQk8yGvfXcvPpVuiNRbl5WwJzJOdxyXhHpqVrrRMTvFOgBdqCpg92H245td/ZE+MJvNjN1XCa//Nh55I/Xw7JEgkSBHkDOOX7x8l6+u66CrnC0374xmak8dMdShblIACnQA+j1fQ18Y205mWkpPPyRZWT0mU6ZmTuKvJxMD6sTkeGiQA+Yju4I//bH3kfVv3zvexivq1REkoYCPUC2HWzmfz62ifKaZi4/c7LCXCTJKNAD5L615ZTXNPO5K+bwr5fO9rocERlhCnSfi0Qdj7y6lwee2UFjew+fv3IOd19W7HVZIuIBBbrPffWJLTz62j4WThvDLecVctsF070uSUQ8okD3sc6eCE9uqqE4bzSPfPx8xmjFIJGkpkD3sfU76mjpCvPDa+YrzEUE3e/tY3/cfIDx2WksnzXR61JEJAEo0H3qzeom/rS5hhsWTSMtpF+jiGjKxTfq27rZUdtybPvBFyoZm5XGPZfP8bAqEUkkCnQf2FnXyk2rXqa+rbtf+2feU8zYLM2di0gvBboP/Ndzb9HVE2HVrUsYk9X7KwuZsahwvMeViUgiUaAnsD9s3M+zFYf446YaPnHRTFYsnOJ1SSKSwBToCaonEuUba8uJRB2Xn5mnuXIROSkFeoL6+1t1NLb3sOrWJRqZi0hcFOgJpqmjhz2H2/jy429SOCGbS+bmel2SiPiEAj2BOOe4adVL7KhtJSM1hcc/tZTMNC3YLCLxUaAniHBsznxHbSt3XjSTfykpYHbeaK/LEhEfUaAngO5wlE/96nWerailcEI2n35PMaMz9KsRkXdGqZEAHiur4tmKWr5x7Xw+/K4ZXpcjIj6lh4B4LBp1/OzF3ZydP5Y7lk/3uhwR8TGN0D30VPlBnnhjPzvr2vjBB87FzLwuSUR8TIHugYNNnXz04VLKa5oZl53G0unjufqsM7wuS0R8ToE+Qjq6I3xj7Rb2HG5n95E22rrCfPqy2XzyktlkpevSRBE5fQr0YRSNOv685SCHWjp5ftsh/rvyMMumT2DB1DF85j3FeriWiAwpBfowqKpv5+5H32BTVeOxtlCK8a3rF3Lr+UUeViYiQaZAH2Ib9jbwgf/7MuGo42MXzuCs/LFcVJxLemoKo3RtuYgMIyXMEPrbjjo++vNSRmWksvq2JZw3U2t9isjIUaAPkR//dSf3/2UbhROyeeTj55E/PtvrkkQkycR1Y5GZrTCz7WZWaWb3DrL/ejPbbGYbzazMzC4c+lIT14a99fznU9tYVDiOX35UYS4i3jjpCN3MQsCDwBVANVBqZmudc1v7HPYcsNY558zsbOAxYN5wFJxo2rrCfPbXm5g6LotffGQZOZla41NEvBHPCH0ZUOmc2+Wc6wbWANf3PcA51+qcc7HNUYAjSTzy6j721bfzwE3nKMxFxFPxBPo0oKrPdnWsrR8zu9HMtgF/Aj4y2BuZ2Z2xKZmyurq6U6k3YUSijsOtXawp3ceiwnE6ASoinosn0Ad7wMhxI3Dn3O+dc/OAG4BvDfZGzrnVzrkS51xJbq6/V+K55aevUPLtZ9lZ18ZdF8/yuhwRkbiucqkGCvps5wM1JzrYObfezGaZ2STn3OHTLTCRhCNRmjp6+NoftvDKrnr+pSSf686ZxoXFk7wuTUQkrkAvBYrNbAawH1gJfLDvAWY2G9gZOym6GEgHjgx1sV5yzvG+H7/EpuomAC6cPYn7rltAdrqu/BSRxHDSNHLOhc3sbuApIAQ85JwrN7O7YvtXAe8DbjezHqAD+ECfk6S+93T5QdaUVrGpuomblxVy7TlnsHyWRuUikljMq9wtKSlxZWVlnnz2O9ETiXLh/c/TFY5y1rSx/OT2Ei3cLCKeMbMNzrmSwfZpvuBtOOf41pNbqW3u4mcfXsqlc/O8LklE5IQU6CdQ09jB9/6yjSc21nDnRTO5ZI6/r8oRkeBToA+itSvM+3/8EjVNndyxfDpfvmqelocTkYSnQB/EA09v50BzJz+5vYQr5k/2uhwRkbgo0PuIRh0/eHYHP3txD7dfUKQwFxFfietpi8nib2/V8V/PV5Kbk8EX3jvX63JERN4RjdBj3qxu4uMPl5GTmcrfv3ipLk0UEd/RCB1obO/mrl9uIJRi3HvVPIW5iPiSRujAA0/v4FBLJ7+9aznnFIzzuhwRkVOS9CP02uZOfvXqXq5aeIbCXER8LakDPRJ1vH/VS0Qd3LBoqtfliIiclqSdcmlo6+ap8oNU1XfwtWvm67Z+EfG9pAz00j313LTqZQDOyR/Lhy4o0p2gIuJ7SRfoVfXt3P3I64zJTOVr18znyvlTSA0l9cyTiARE0gX6A09vp7a5i3suL+amkoKTf4OIiE8k1dC0tSvMn7cc5NbzC7nn8jlelyMiMqSSKtBLd9fTFY5y1cIzvC5FRGTIJVWgv1h5mPRQCkuKxntdiojIkEuaQK9v6+axsireXTxJt/aLSCAlTaC/sO0QzZ1hPnN5sdeliIgMi6QJ9B21LaSHUph/xhivSxERGRZJFegzc0fpmnMRCaykSLeDTZ1sqm6ieHKO16WIiAybpAj0H/+1krauMB9513SvSxERGTZJEehlextYOn0Ciwp1uaKIBFfgA72tK0zFgWYW69pzEQm4wAf6xqpGog7dTCQigRf4QN+wtwEzOFerEYlIwAU60J1zrN9Rx5y8HMZmpXldjojIsAp0oL+88whlextYuUyPyRWR4At0oD9VfpCstBA3Lyv0uhQRkWEX6EBf/9Zhzp85QQ/jEpGkENhAr23uZPfhNt41e5LXpYiIjIjABnrpnnoAlk6f4HElIiIjI7iBvrue7PQQC6bq6YoikhyCG+h7GlhUOE5PVxSRpBHItGvu7KHiYDMlRZpuEZHkEVegm9kKM9tuZpVmdu8g+28xs82xr5fM7JyhLzV+r+9twDlYNkOBLiLJ46SBbmYh4EHgKmA+cLOZzR9w2G7gYufc2cC3gNVDXeg7UbqnnlCK6XZ/EUkq8YzQlwGVzrldzrluYA1wfd8DnHMvOecaYpuvAPlDW+Y78+quehZMHcOojFQvyxARGVHxBPo0oKrPdnWs7UQ+Cvx5sB1mdqeZlZlZWV1dXfxVvgO1zZ1s2NfApXPzhuX9RUQSVTyBboO0uUEPNLuU3kD/0mD7nXOrnXMlzrmS3Nzc+KuMk3OOf19XgXNw7TlTh/z9RUQSWTxzEtVA36db5QM1Aw8ys7OBnwJXOeeODE1578zGqkae2FjD3ZfOZnbeaC9KEBHxTDwj9FKg2MxmmFk6sBJY2/cAMysEHgduc87tGPoy4/NYWRWj0kN84uKZXpUgIuKZk47QnXNhM7sbeAoIAQ8558rN7K7Y/lXA14GJwI/MDCDsnCsZvrIH92LlEZbPnkROpp59LiLJJ67LQJxz64B1A9pW9Xn9MeBjQ1vaO7O/sYN99e3cfkGRl2WIiHgmMHeKfnddBWkh4/IzJ3tdioiIJwIR6M9ureVPmw9wz+VzmD5plNfliIh4IhCB/sTG/eTmZPCJi3QyVESSl+8DPRJ1/P2tw1w8J1dPVhSRpOb7BNxV10pTRw8XzJzodSkiIp7yfaBvPdAMwIJpWshCRJJbIAI9PZTCrFzdGSoiyc33gb5xXyNzp+SQpvlzEUlyvk7Bju4Ib+xr5IJZmj8XEfF1oJftrac7ElWgi4jg80B/aecRUlOMZdO11JyIiO8D/dyCcVqZSEQEHwe6c47tB5s5O1/rhoqIgI8D/UhbN509UQomZHldiohIQvBtoFc3dACQPz7b40pERBKDjwO9HYD88Rqhi4iArwO9d4Q+TYEuIgL4ONCPtHaRmZbCGC03JyIC+DjQG9t7GJeV7nUZIiIJw7eB3tTRw9gsjc5FRI7ybaA3dvQwNluBLiJylG8DvVkjdBGRfnwb6E0dPYxToIuIHOPbQG9s72GcplxERI7xZaB3hSN09EQ05SIi0ocvA72powdAgS4i0ocvA721MwzAGAW6iMgx/gz0rt5AH5Wu56CLiBzlz0CPjdBHZyrQRUSO8megx0boo7VSkYjIMQp0EZGA8GWgtx2dQ1egi4gc48tAb4kFeo7m0EVEjvFloLd1hQmlGBmpvixfRGRY+DIRWzvDjM5Ixcy8LkVEJGH4MtBbusI6ISoiMoAvA71NgS4icpy4At3MVpjZdjOrNLN7B9k/z8xeNrMuM/v80JfZX1tXhOyM0HB/jIiIr5x0mGtmIeBB4AqgGig1s7XOua19DqsHPg3cMCxVDtDeHdZt/yIiA8QzQl8GVDrndjnnuoE1wPV9D3DOHXLOlQI9w1Djcdq7I2SmaYQuItJXPIE+Dajqs10da/NMR0+E7HQFuohIX/EE+mDXBrpT+TAzu9PMysysrK6u7lTeAoCObgW6iMhA8QR6NVDQZzsfqDmVD3POrXbOlTjnSnJzc0/lLYDeQM9SoIuI9BNPoJcCxWY2w8zSgZXA2uEt68Scc7T3RMjSHLqISD8nvVTEORc2s7uBp4AQ8JBzrtzM7ortX2VmU4AyYAwQNbN7gPnOueahLrg7EiUSdZpyEREZIK5r/5xz64B1A9pW9Xl9kN6pmGHX2R0FIEuXLYqI9OO7O0Xbe3qftKgRuohIf/4L9O4IgObQRUQG8F2gdxwNdI3QRUT68V+g9/QGuqZcRET6812gH51yUaCLiPTnu0Dv6O49KapnuYiI9Oe7QM/NyeDqs6YwcVSG16WIiCQU313MvaRoAkuKJnhdhohIwvHdCF1ERAanQBcRCQgFuohIQCjQRUQCQoEuIhIQCnQRkYBQoIuIBIQCXUQkIMy5U1rv+fQ/2KwO2HuK3z4JODyE5fiB+pwc1OfkcDp9LnLODboos2eBfjrMrMw5V+J1HSNJfU4O6nNyGK4+a8pFRCQgFOgiIgHh10Bf7XUBHlCfk4P6nByGpc++nEMXEZHj+XWELiIiAyjQRUQCwneBbmYrzGy7mVWa2b1e1zNUzOwhMztkZlv6tE0ws2fM7K3Yn+P77Pty7Gew3cze603Vp8fMCszsBTOrMLNyM/tMrD2w/TazTDN7zcw2xfr8zVh7YPsMYGYhM3vDzJ6MbQe6vwBmtsfM3jSzjWZWFmsb3n4753zzBYSAncBMIB3YBMz3uq4h6ttFwGJgS5+27wH3xl7fC9wfez0/1vcMYEbsZxLyug+n0OczgMWx1znAjljfAttvwIDRsddpwKvA+UHuc6wfnwMeAZ6MbQe6v7G+7AEmDWgb1n77bYS+DKh0zu1yznUDa4DrPa5pSDjn1gP1A5qvBx6OvX4YuKFP+xrnXJdzbjdQSe/Pxleccwecc6/HXrcAFcA0Atxv16s1tpkW+3IEuM9mlg/8E/DTPs2B7e9JDGu//Rbo04CqPtvVsbagmuycOwC94QfkxdoD93Mws+nAInpHrIHud2z6YSNwCHjGORf0Pv8A+CIQ7dMW5P4e5YCnzWyDmd0ZaxvWfvttkWgbpC0Zr7sM1M/BzEYDvwPucc41mw3Wvd5DB2nzXb+dcxHgXDMbB/zezBa+zeG+7rOZXQMccs5tMLNL4vmWQdp8098B3uWcqzGzPOAZM9v2NscOSb/9NkKvBgr6bOcDNR7VMhJqzewMgNifh2Ltgfk5mFkavWH+K+fc47HmwPcbwDnXCPwVWEFw+/wu4Doz20PvFOllZvZLgtvfY5xzNbE/DwG/p3cKZVj77bdALwWKzWyGmaUDK4G1Htc0nNYCH4q9/hDwhz7tK80sw8xmAMXAax7Ud1qsdyj+/4AK59z3++wKbL/NLDc2MsfMsoDLgW0EtM/OuS875/Kdc9Pp/fv6vHPuVgLa36PMbJSZ5Rx9DVwJbGG4++31meBTOHN8Nb1XQ+wEvuJ1PUPYr0eBA0APvf9afxSYCDwHvBX7c0Kf478S+xlsB67yuv5T7POF9P5v5WZgY+zr6iD3GzgbeCPW5y3A12Ptge1zn35cwj+ucgl0f+m9Em9T7Kv8aFYNd79167+ISED4bcpFREROQIEuIhIQCnQRkYBQoIuIBIQCXUQkIBToIiIBoUAXEQmI/w/CbtzryaXpOAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# hard constraints\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import cvxpy as cp\n",
    "\n",
    "N = 1\n",
    "T = 500\n",
    "m = 20  # number of devices\n",
    "n = 5  # number of servers\n",
    "y_max = 110\n",
    "y_min = 90\n",
    "rs = 0.1\n",
    "beta = np.ones(n)\n",
    "beta[0] = 3\n",
    "beta[1] = 3\n",
    "seed = 0\n",
    "np.random.seed(seed)\n",
    "\n",
    "def optimization(m, n, wS, y, mu, BS, beta, CD):\n",
    "    x = cp.Variable((n + 1, m), nonneg=True)\n",
    "    obj = beta @ cp.inv_pos(cp.sqrt(x[1:, :] @ y + wS))\n",
    "\n",
    "    constraints = [0 <= x, x <= 1,\n",
    "                   x[1:, :] @ y <= BS,\n",
    "                   cp.sum(x, 0) == 1]\n",
    "    for i in range(m):\n",
    "        constraints.append((y[i] * mu[i, :] @ x[1:, i]) <= CD[i])\n",
    "\n",
    "    prob = cp.Problem(cp.Minimize(obj), constraints)\n",
    "    prob.solve()  # Returns the optimal value.\n",
    "    return x, prob\n",
    "\n",
    "def oracle(y, mu):\n",
    "    x, prob = optimization(m, n, wS, y, mu, BS, beta, CD)\n",
    "    return x.value, prob.value, prob.status\n",
    "\n",
    "def f(x, y):\n",
    "    return beta.dot(1/np.sqrt(x[1:, :].dot(y) + wS))\n",
    "\n",
    "def f_drop(x, y, mu, CD, BS):\n",
    "    yp = np.minimum(y, CD/ np.diag(mu @ x[1:, :]))\n",
    "    dp = np.minimum(x[1:, :].dot(yp), BS)\n",
    "    return beta.dot(1/np.sqrt(dp + wS))\n",
    "\n",
    "reg = np.zeros((N, T))\n",
    "\n",
    "#records all y, x_opt, x_t #yuhang yao\n",
    "y_N_T = np.zeros((N, T, m)) #yuhang yao\n",
    "x_opt_N_T = np.zeros((N, T, n + 1, m)) #yuhang yao\n",
    "x_t_N_T = np.zeros((N, T, n + 1, m)) #yuhang yao\n",
    "j_N_T = np.zeros((N, T, m))#yuhang yao\n",
    "BS_N = np.zeros((N, n))#yuhang yao\n",
    "\n",
    "for u in range(N):\n",
    "    wS = np.random.randint(15, 25, n)\n",
    "    \n",
    "#     BS = np.random.uniform(y_min*5, y_max*5, n)\n",
    "#     CD = np.random.uniform(y_min/2, y_max/2, m)\n",
    "    BS = np.random.uniform(y_min*10, y_max*10, n)\n",
    "    BS_N[u] = BS\n",
    "    CD = np.random.uniform(y_min*0.7, y_max*0.7, m)\n",
    "    mu = np.random.rand(m, n)\n",
    "    mu[:,0] = 0.5\n",
    "    mu[:,1] = 0.5\n",
    "    mu[:,2] = 0.8\n",
    "    mu[:,3] = 0.8\n",
    "    mu[:,4] = 0.8\n",
    "    # trace_gen = Trace(m, n, seed + u)\n",
    "    # mu = trace_gen.avg()\n",
    "    # mu = np.random.rand(m, n)\n",
    "    # mu_hat = np.zeros_like(mu)  # empirical mean\n",
    "    mu_hat = np.ones_like(mu)\n",
    "    T_ij = np.ones_like(mu)  # total number of times arm (i,j) is played\n",
    "    for t in range(T):\n",
    "        y = np.random.uniform(y_min, y_max, m).astype(int)\n",
    "        x_opt, f_opt, status = oracle(y, mu)\n",
    "        if 'optimal' not in status:\n",
    "            print('Solution infeasible 1')\n",
    "            break\n",
    "\n",
    "        rho_ij = np.sqrt(3 * np.log(t + 1) / (2 * T_ij)) * rs\n",
    "        mu_bar = np.maximum(mu_hat - rho_ij, 0) # LCB\n",
    "        x_tmp, f_tmp, status = oracle(y, mu_bar)\n",
    "        if 'optimal' not in status:\n",
    "            print('Solution infeasible 2')\n",
    "            break\n",
    "        \n",
    "        # mapping from x_tmp to x_t (new)\n",
    "        x_t = np.zeros((n + 1, m))\n",
    "        \n",
    "        for i in range(m):\n",
    "            cost = y[i] * mu[i,:].dot(x_tmp[1:, i])\n",
    "            if cost > CD[i]:\n",
    "                x_t[1:, i] = mu_bar[i] / mu[i] * x_tmp[1:, i]\n",
    "                x_t[0, i] = 1-np.sum(x_t[1:, i])\n",
    "                if x_t[0, i] < 0:\n",
    "                    x_t[0, i] = 0\n",
    "                    x_t[1:, i] = x_t[1:, i]/np.sum(x_t[1:, i])\n",
    "            else:\n",
    "                x_t[:, i] = x_tmp[:, i]\n",
    "\n",
    "        f_t = f(x_t, y)\n",
    "        \n",
    "        # sample j based on x_t[i], observe c_ij, update mu_hat[i,j]\n",
    "        # c = trace_gen.generate()\n",
    "        for i in range(m):\n",
    "            j = np.random.choice(n+1, p=x_t[:, i])\n",
    "            j_N_T[u, t, i] = j #yuhang yao\n",
    "            if j != 0:\n",
    "                j -= 1\n",
    "                c_ij = int(np.random.rand() < mu[i, j])\n",
    "                # a = np.random.rand() * 3\n",
    "                # c_ij = np.random.beta(a, a * (1-mu[i, j])/mu[i, j]) # beta distribution\n",
    "                # c_ij = c[i, j]  # trace\n",
    "                T_ij[i, j] += 1\n",
    "                mu_hat[i, j] += (c_ij - mu_hat[i, j]) / T_ij[i, j]\n",
    "\n",
    "        # calculate regert\n",
    "        reg[u, t] = f_t - f_opt\n",
    "        y_N_T[u, t] = y#yuhang yao\n",
    "        x_opt_N_T[u, t] = x_opt#yuhang yao\n",
    "        x_t_N_T[u, t] = x_t#yuhang yao\n",
    "        \n",
    "plt.plot(np.cumsum(reg, axis=1).T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([605.70735018, 605.70937788, 291.19441505, 291.19444473,\n",
       "       291.19441143])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_t[1:, :].dot(y) + wS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([605.71613982, 605.70536136, 291.19407864, 291.1941644 ,\n",
       "       291.19024448])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_opt[1:, :].dot(y) + wS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1240.6749809965972"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum((y*x_t[1:]).T*mu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c3a624f280>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3yU133n8c9PAnEXkpDAoAsgIzCYm7G4OL7E1xjsNDTdboodx7GThrg1uW03WSfpNk376m7aV9uNs3VNiE0cb9I4qRO3NKG2c7FD4hiMuIO5CRmQkEASQiMJIXSZ3/4xD2QsC2uMJY30zPf9es1LM89zRvM7Ar48OnPmHHN3REQkvNKSXYCIiPQvBb2ISMgp6EVEQk5BLyIScgp6EZGQG5bsAnqSm5vr06ZNS3YZIiJDxrZt2+rdPa+nc4My6KdNm0ZZWVmyyxARGTLM7NilzmnoRkQk5BT0IiIhp6AXEQk5Bb2ISMgp6EVEQk5BLyIScgp6EZGQU9CLiCRZw9l2/n3nCR5/+Ui/fP+EPjBlZsuBR4F04Al3/1q38xacvwtoBR5w9+3Buc8Bfww4sAd40N3b+qwHIiJDTDTq7DkR4eWDdbx0sJZdVY24w+TxI/nEjdMZlt631+C9Br2ZpQOPAXcAVcBWM9vg7q/HNVsBlAS3pcDjwFIzywc+Dcxx93Nm9kNgFfBUn/ZCRGQQcnfqW9qpOtNK5ZlzVDa08np1E5srTnP6bDtmML8gi8/cVsLNsyYyL3886WnW53UkckW/BCh39woAM3sGWAnEB/1K4GmPbVe12cyyzGxy3GuMMrMOYDRQ3WfVi4gkWXNbB5UN5zje0BoL9IbfhXrVmXOc6+h6U/uC7FHcUJLLzbPyuKkkjwljR/R7jYkEfT5QGfe4ithVe29t8t29zMz+HjgOnANedPcXe3oRM1sNrAYoKipKrHoRkX7WFXWqG2PBfTy4HWtovfi4sbXjTe3HjRhGQc5opueO4aaZeRRmj6IgezSFOaMpyB7FmBEDv8RYIq/Y0+8R3Tea7bGNmWUTu9qfDjQC/2pm97n7d9/S2H0dsA6gtLRUG9mKyICKnOvgQE0TB042c+hUM8eDMK86c47O6O8iaViaUZA9isKc0dw1bzJFOaMpzB4d+5ozivGjhhN723LwSCToq4DCuMcFvHX45VJtbgfecPc6ADP7MfAe4C1BLyIyUBpb29l7ook9JyLsPRFhb3WEY6dbL54fP2o4UyeMZm7++IthHgvy0UweP7LP3yztb4kE/VagxMymAyeIvZl6b7c2G4A1wfj9UiDi7jVmdhxYZmajiQ3d3AZo/WERGTBNbR3sqYqwuyrCnhON7DkRobLh3MXzBdmjmJc/ng+VFnL1lExmT85k4rgRg+6q/N3oNejdvdPM1gAvEJteud7d95nZQ8H5tcBGYlMry4lNr3wwOLfFzJ4FtgOdwA6C4RkRkb52rr2LfdURdlVF2FPVyO6qCBX1Zy+eL8oZzfz8LO5dMpV5+eOZm59J1uiMJFY8MCw2UWZwKS0tdW08IiJvp70zysGTzew+0cjuygi7qho5XNtCVzCefkXmSOYVjGdBwXjmF2Qxv2B8qEPdzLa5e2lP5wblDlMiIvG6ok5FXQu7qiLsrmpkV1WE/TVNtHdGAcgaPZz5BVncMWfSxVCflDkyyVUPHgp6ERlU3J2aSBs7KxvZVdnIzspG9p6IcLY9Nh99TEY6c/PH88B7pjG/YDzz87MozBkVqjH1vqagF5GkauvoYu+JCNuOnWH78TPsON5IbfN5ADLS05gzJZM/vLaA+QVZLCgcT3HuWNL64dOjYaagF5EBVRM5x/ZjjReDfV91hI6u2Lj61AmjuX5GLgsLs1hYmMXsyZlkDBtaUxkHIwW9iPSbjq4o+2ua2HbsTCzYj52hOhJb03DEsDQWFGTx8RuKWVSUxaKp2eQOwHIAqUhBLyJ9prW9k23HzrClooGyYw3sqoxcXOtlyviRLJqazSemZrOoKFtX6wNIQS8il+3s+U7Kjp1hS8VpNlecZndVhM6ok55mXD0lk1VLCrl2ajbXTs1m8vhRyS43ZSnoRSRhZ893svVoA5srGthccZo9JyJ0RZ1hacb8gvGsvqmYZcUTuHZqdlIW75Ke6U9CRC6praOL7cfOsLniNK8cOc2uykY6o87wdGNBQRZ/8t4rWVqcw7VTsxmdoTgZrPQnIyIXnWvvYvvxWLBvqWhgZ2Uj7V1R0oINMlbfVMx1V06gdGoOozLSk12uJEhBL5LCWs53Una0gS1vNLAlGIrp6HLSDOblj+fB66exrHgCpdOyGTdyeLLLlcukoBdJIc1tHZQFQzGbKxrYGzfGPq9gPB+/oZil03MU7CGjoBcJsQtvnr7aLdiHpxsLCzXGnir0JysSIm0dXew43sirR+r57ZHT7Ix78/SawmwevvlKlhVP4JqibI2xpxAFvcgQ5u4cPNXMpkN1bDpUz9ajDZzv1Jun8mYKepEhpqmtg1cO1/PywTp+daiOk02xJQVmThrLh5dO5foZE1g8PYdMjbFLIKGgN7PlwKPEdph6wt2/1u28BefvIrbD1APuvt3MZgE/iGtaDPyFu3+9L4oXSQVtHV3sq25i69EGXj5YS9nRM3RGnXEjh3FjSS7vnZnHTTPz9MlTuaReg97M0oHHgDuIbQK+1cw2uPvrcc1WACXBbSnwOLDU3Q8CC+O+zwnguT7tgUiIuDsnGs9dXABs2/EzHKhppjPYNWn25ExW31TMzbMmsqgoa8htUi3JkcgV/RKg3N0rAIINwFcC8UG/EnjaY/sSbjazLDOb7O41cW1uA464+7E+ql1kyGvvjLKv+ndrsW87doZTTbG12EdnpLOwMDbOvrAwi4VFWUwcp12T5J1LJOjzgcq4x1XErtp7a5MPxAf9KuD7l3oRM1sNrAYoKipKoCyRoed8Zxfbjp5h0+F6th1rYHdVhPPBdniFOaO4LlgnZtHUbGZNGqcrdukTiQR9T1u5dN9R/G3bmFkG8AHgi5d6EXdfB6yD2ObgCdQlMui5O0fqzvKbw3X8pjw25bG1vYvh6cbc/PF8ZNnUi6s7TtQep9JPEgn6KqAw7nEBUP0O26wAtrv7qcspUmQoaW3v5JXy0/zyQC2/Olh7caONqRNG8weL8rl55kSuu3KCVneUAZPI37StQImZTSf2Zuoq4N5ubTYAa4Lx+6VApNv4/D28zbCNyFB3vrOLTYfq2bCrmp+/fopzHV2MHTGMG2bksubWPG4syaUwZ3Syy5QU1WvQu3unma0BXiA2vXK9u+8zs4eC82uBjcSmVpYTm1754IXnm9loYjN2Ptn35YskT0dXlFfK6/np7hpe2HeSprZOskcP54OL8rl73mQWT8vRDkoyKCT0u6O7byQW5vHH1sbdd+DhSzy3FZjwLmoUGTQ6u6K8WnGan+6u4fl9J2ls7WDciGHcMWcSv7dwCjfMyGW43kCVQUaDhCK96OyKsuWNBn4SXLk3nG1nbBDud8+bzI0zcxkxTMsLyOCloBfpQVfUee2NBn66p5rn956kvqWd0Rnp3D57EnfPn8x7Z+YxcrjCXYYGBb1IIBp1th5t4Kd7ati45yT1LecZNTyd22ZP5P3zJ3PzrIkKdxmSFPSS0s53drH1jTP8fP8pNu6pobb5PCOHp3HrVRO5e94UbrkqT+u0y5Cnv8GSciKtHfziwCle3HeKTYfraG3vYsSwNG6elcf750/h1qsmao67hIr+NkvouTvltS28Ul7Pz/fXsrniNJ1RZ1LmCD54TT63XjWR91yZq/XaJbQU9BJKdc3n+fXhOjYdquM35aepb4ktFFacO4Y/vrGYO6+exIKCLNLSelq9QyRcFPQSCm0dXWw/foZNh+rZdKiO12uaAMgZk8ENM3K5fsYErivOpTBnFLHtE0RSh4JehpzOriiHTrWwu6qRXVURdlU2cuhUbM32YWnGoqnZfP7OWdxUksfVUzJ11S4pT0Evg1514zm2Hm1gV2WE3VWN7K2O0NYRW9o3c+Qw5hdk8cn3FrOwMJtlxTmM0xZ6Im+ioJdBp7G1nVePnOaVI/W8Un6aN+rPAjByeBpXTxnPvUumsqBwPPMLspg2YbSGYkR6oaCXpOuKOjuOn+EXB2p5pbyePSciuMOYjHSWFU/gvmVTWVaco404RC6Tgl6S4tjps7xSfprfHqnn1SOnOX22PTa+XpTNZ2+byfUzJrCgMEsLhIn0AQW9DJjy2mZ+vP0EL75+ivLaFgAmZY7gxpJcbp09iVtm5Wl8XaQfKOil3+2sbOSfXyrnxddPkZ5mLCvO4b6lRdw4M4/i3DEaYxfpZwp66Rfuzm+PnOafXy7nlfLTjB81nE/fVsJHr5vKhLEjkl2eSEpJKOjNbDnwKLEdpp5w9691O2/B+buI7TD1gLtvD85lAU8Ac4ltGP4xd3+1z3ogg0o06vxs/yn++eUj7KpsZOK4EXz5rtncs7SIsVo/RiQpev2XZ2bpwGPEtgOsAraa2QZ3fz2u2QqgJLgtBR4PvkLsP4Dn3f0PzSwD0MaZIdTRFeU/dlXz+MtHOFzbQlHOaP7XB+fxB4vytbSvSJIlcom1BCh39wqAYAPwlUB80K8Eng62FNxsZllmNhk4C9wEPADg7u1Ae9+VL8nW1tHFD8sq+eavKjjReI6rrhjHo6sWcve8yZoKKTJIJBL0+UBl3OMqfne1/nZt8oFOoA74tpktALYBn3H3s91fxMxWA6sBioqKEq1fkqSjK8p3fnuUtb86Qn1LO9dOzeavVl7NrVdN1JurIoNMIkHf079aT7DNMGAR8Cl332JmjwKPAP/zLY3d1wHrAEpLS7t/fxkk3J2XD9XxtY0HOHiqmRtLcllzywyWTM9RwIsMUokEfRVQGPe4AKhOsI0DVe6+JTj+LLGglyFoZ2UjX/vP/WyuaKAoZzTrPnIt77v6imSXJSK9SCTotwIlZjYdOAGsAu7t1mYDsCYYv18KRNy9BsDMKs1slrsfBG7jzWP7MgRU1LXw9y8eZOOek0wYk8FXP3A19ywpImOYxuBFhoJeg97dO81sDfACsemV6919n5k9FJxfC2wkNrWynNj0ygfjvsWngO8FM24qup2TQexUUxtf//lhflhWyYhhaXz6thJW31SsaZIiQ4zFJsoMLqWlpV5WVpbsMlLWyUgb3/p1Bd/bcoyuqHPvkiLW3FpC3jh90ElksDKzbe5e2tM5XZrJRZUNrTz2Ujk/2l5F1GHlgil85vYSpk4Yk+zSRORdUNALtU1t/N9flvPM1uOkmbFqcRGrbyqmMEefbRMJAwV9Cjtztp21m47wnd8epbPLWbWkkE/dWsKkzJHJLk1E+pCCPgWdPd/J+t+8wbpNFbS0d/LBhfl89vaZFE3QFbxIGCnoU8zPXj/FV/59L9WRNt43ZxJ/9r5ZzLpiXLLLEpF+pKBPEScjbfzlhn08v+8ksyaN49l7rqF0Wk6yyxKRAaCgD7lo1PnulmP83fMH6eiK8oXls/jEjcXaok8khSjoQ+xUUxv//V938evD9dxYksvf/P48jcOLpCAFfUj9dHcNX3puD+2dUf769+dy39IiLTomkqIU9CHT1tHFl368hx/vOMGCgvH8nz9aSHHe2GSXJSJJpKAPkea2Dj7xdBmbKxr49K0z+NRtJRqLFxEFfVhEWjt44KnX2FMV4dFVC1m5MD/ZJYnIIKGgH+LcnQ27qvnrn+wncq6dxz68iDu1RryIxFHQD2H1Lef58nN7eGHfKRYUZvHtBxYzr2B8sssSkUFGQT9EPb+3hi89t5eWtk4eWXEVn7ixmPQ0zaoRkbdKKOjNbDnwKLGNR55w9691O2/B+buIbTzygLtvD84dBZqBLqDzUuslS2IirR18ZcNe/m1nNXPzM/nHDy1k5iQtYSAil9Zr0JtZOvAYcAexvWG3mtkGd4/fEnAFUBLclgKPB18vuMXd6/us6hT10sFaHvnRbk63tPO522fyp7dcqVk1ItKrRK7olwDl7l4BEOwLu5I37/26EnjaY9tVbTazLDObfGHfWHl33J1/ePEQ//RSOTMnjeXJjy5mbr7G4kUkMYlcDuYDlXGPq4JjibZx4EUz22Zmqy+30FQVjTpf2bCPf3qpnFWLC/mPT92gkBeRdySRK/qe3uHrvtHs27W53t2rzWwi8DMzO+Dum97yIrH/BFYDFBUVJVBW+HV2Rfn8s7t5bscJPnlTMY+suErLGIjIO5bIFX0VUBj3uACoTrSNu1/4Wgs8R2wo6C3cfZ27l7p7aV5eXmLVh1h9y3k++u3XeG7HCT5/5yyFvIhctkSCfitQYmbTzSwDWAVs6NZmA3C/xSwDIu5eY2ZjzGwcgJmNAd4H7O3D+kPptTcauOvRX1N29Ax/94fzefiWGQp5EblsvQ7duHunma0BXiA2vXK9u+8zs4eC82uBjcSmVpYTm175YPD0ScBzQUgNA/7F3Z/v816EREdXlMdfPsKjvzhMYfYonnpwCXOmZCa7LBEZ4iw2UWZwKS0t9bKysmSXMaDKa1v49Pd38HpNE7+3YAp/88G5ZI4cnuyyRGSIMLNtl/qckj4ZOwi8Ul7Pn3x3G8PT01h737Usn6u1akSk7yjok+yZ147z5/+2l+K8MTz50cUU5mgHKBHpWwr6JIlGnb99/gDf3FTBjSW5PPbhRRqqEZF+oaBPgtb2Tj73g528sO8UH15axFc/cDXDtJSBiPQTBf0Ai7R2cP+3X2N3VSN/fvdsPn7DdE2dFJF+paAfQI2t7dz35BYOnWxh7X3XaoMQERkQCvoB0nC2nfue2EJ5XQvf/Mi13HLVxGSXJCIpQkE/AE63nOfDT2yhov4s37q/lPfO1BIPIjJwFPT9rL7lPB/+1haOnj7L+o8u5oaS3GSXJCIpRkHfj9o7o3zsqa0cazjLtx9YzHtmKORFZOAp6PvRN35xmN1VEdbed61CXkSSRpO3+8m+6ghrf3WEP1iUryUNRCSpFPT94MzZdj77zE6yRmfwF++fk+xyRCTFaeimjx0+1cyDT22ltuk8Tz24mKzRGckuSURSnIK+D+2rjvCRJ18jPc344UPXsbAwK9kliYgo6PvK7qpGPvLka4zJSOd7n1jG9NwxyS5JRARIcIzezJab2UEzKzezR3o4b2b2jeD8bjNb1O18upntMLOf9FXhg8meqggf/tYWMkcN4wefvE4hLyKDSq9Bb2bpwGPACmAOcI+ZdX+HcQVQEtxWA493O/8ZYP+7rnYQqm1q46HvbiNz1HB+sPo6rScvIoNOIlf0S4Byd69w93bgGWBltzYrgac9ZjOQZWaTAcysALgbeKIP6x4U6lvOc+8TW2hsbefx+xYxJWtUsksSEXmLRII+H6iMe1wVHEu0zdeBLwDRt3sRM1ttZmVmVlZXV5dAWclV33Kee7+1maozrTz5wGLmF+iNVxEZnBIJ+p4WS+++o3iPbczs/UCtu2/r7UXcfZ27l7p7aV7e4F70q7W9kwe/vZXjDa2sf2Axy4onJLskEZFLSiToq4DCuMcFQHWCba4HPmBmR4kN+dxqZt+97GoHgY6uKJ/6lx3sq47wT/cs4j1XamkDERncEgn6rUCJmU03swxgFbChW5sNwP3B7JtlQMTda9z9i+5e4O7Tguf90t3v68sODKRo1PnCs7v5xYFa/mrlXG6fMynZJYmI9KrXefTu3mlma4AXgHRgvbvvM7OHgvNrgY3AXUA50Ao82H8lJ4e781c/eZ3ndpzg83fO4r5lU5NdkohIQhL6wJS7byQW5vHH1sbdd+DhXr7Hy8DL77jCQeJbv67gqd8e5WPXT+dPb74y2eWIiCRMi5ol4JcHTvG///MAd8+bzJ/fPVubeYvIkKKg78WRuhY+8/2dzJmcyd//1wWkpSnkRWRoUdC/jaa2Dj7xdBkZw9JYd38pozLSk12SiMg7pkXNLqEr6nz2mZ0cP93K9/54Kfn61KuIDFG6or+E//fqUX55oJavfOBqluoDUSIyhCnoexA518E3N1WwZHoO9y0tSnY5IiLvioK+m9b2Tj721FbqW87z+TtnaYaNiAx5Cvo47Z1R/uS729lx/AyPrrqGxdNykl2SiMi7pjdjA9Go87kf7uRXh+r42/8yj7vmTU52SSIifUJX9IHvvXacn+6u4X8sv4o/WqxxeREJDwU90NzWwT+8eJDrZ0zgofcWJ7scEZE+paAH1v/mKI2tHTyyXMsbiEj4pHzQN7a288SvK7jz6knMKxif7HJERPpcygf9+leO0tLeyefumJnsUkRE+kVKB31X1PnXskpuKsnjqisyk12OiEi/SOmg/015PTWRNv5ocWHvjUVEhqiEgt7MlpvZQTMrN7NHejhvZvaN4PxuM1sUHB9pZq+Z2S4z22dmX+3rDrwbz26rInv0cG6bPTHZpYiI9Jteg97M0oHHgBXAHOAeM5vTrdkKoCS4rQYeD46fB2519wXAQmB5sKds0nVFnU2H6rhjziRGDNPywyISXolc0S8Byt29wt3bgWeAld3arASe9pjNQJaZTQ4etwRthgc376vi3439NU1EznVw/YzcZJciItKvEgn6fKAy7nFVcCyhNmaWbmY7gVrgZ+6+pacXMbPVZlZmZmV1dXWJ1n/Z/mN3NelpxnVXagliEQm3RIK+p08Qdb8qv2Qbd+9y94VAAbDEzOb29CLuvs7dS929NC8vL4GyLl9reyff33Kc5VdfwcRxI/v1tUREki2RoK8C4qelFADV77SNuzcCLwPL33GVfexH20/Q1NbJx26YluxSRET6XSJBvxUoMbPpZpYBrAI2dGuzAbg/mH2zDIi4e42Z5ZlZFoCZjQJuBw70Yf2X5Ufbqpg9OZNFRdnJLkVEpN/1ukyxu3ea2RrgBSAdWO/u+8zsoeD8WmAjcBdQDrQCDwZPnwx8J5i5kwb80N1/0vfdSNzx063srGzkkRVXaV0bEUkJCa1H7+4biYV5/LG1cfcdeLiH5+0GrnmXNfapZ7dVYgYfWDAl2aWIiAyIlPpkbFfUeXZbFTeW5DEla1SyyxERGRApFfSvHjlNdaSND5UWJLsUEZEBk1JB/8qReoanG7fPnpTsUkREBkxKBf3eExFmThrHyOFa8kBEUkfKBL27s/dEhLlTtLmIiKSWlAn6E43nONPawVztIiUiKSZlgn7viSYA5uUr6EUktaRQ0EdITzOuumJcsksRERlQKRP0ZccauOoKvRErIqknJYK+tb2T7ccatfa8iKSklAj6Hccbae+K8h6tPS8iKSglgn7PiQgACwuzklyJiMjAS4mg33siQkH2KLJGZyS7FBGRAZcSQf96TRNXT8lMdhkiIkkR+qCPRp2qhnNMyx2T7FJERJIi9EF/qrmN9q4ohdmjk12KiEhSJBT0ZrbczA6aWbmZPdLDeTOzbwTnd5vZouB4oZm9ZGb7zWyfmX2mrzvQm8qGcwAU5ijoRSQ19Rr0wTaAjwErgDnAPWY2p1uzFUBJcFsNPB4c7wT+zN1nA8uAh3t4br+qbGgFoDBbG42ISGpK5Ip+CVDu7hXu3g48A6zs1mYl8LTHbAayzGyyu9e4+3YAd28G9gP5fVh/r47UtTAszchX0ItIikok6POByrjHVbw1rHttY2bTiO0fu6WnFzGz1WZWZmZldXV1CZSVmJ2VjcyenMmIYVr6QERSUyJBbz0c83fSxszGAj8CPuvuTT29iLuvc/dSdy/Ny8tLoKzedUWdXZWNXFOkD0qJSOpKJOirgMK4xwVAdaJtzGw4sZD/nrv/+PJLfefeqD/L2fYu5hco6EUkdSUS9FuBEjObbmYZwCpgQ7c2G4D7g9k3y4CIu9eYmQFPAvvd/R/7tPIEHDgZ++VBSxOLSCob1lsDd+80szXAC0A6sN7d95nZQ8H5tcBG4C6gHGgFHgyefj3wEWCPme0Mjn3J3Tf2bTd6dvBkM+lpxoyJYwfi5UREBqVegx4gCOaN3Y6tjbvvwMM9PO839Dx+PyAOnGxmeu4YrUEvIikt1J+MPXCyiVkathGRFBfaoG8530llwzlmK+hFJMWFNugPnWoGYNYVWrVSRFJbaIP+wtIH07VqpYikuNAG/clIGwBXjB+Z5EpERJIrtEF/quk8YzLSGTsioYlFIiKhFd6gb25jUqau5kVEwhv0kTYmZo5IdhkiIkkX3qBvbuMKXdGLiIQ36M+c7SB7TEayyxARSbpQBn1X1Gk538n4UcOTXYqISNKFMuhb2joBGDdSQS8iEsqgb2rrACBzpKZWioiEMugj54Kg19CNiEg4g/7CFf04XdGLiIQz6JuDMfpMjdGLiCQW9Ga23MwOmlm5mT3Sw3kzs28E53eb2aK4c+vNrNbM9vZl4W+nKRi60awbEZEEgt7M0oHHgBXAHOAeM5vTrdkKoCS4rQYejzv3FLC8L4pNVJOu6EVELkrkin4JUO7uFe7eDjwDrOzWZiXwtMdsBrLMbDKAu28CGvqy6N5cuKIfqzF6EZGEgj4fqIx7XBUce6dt3paZrTazMjMrq6ureydPfYvmtk7GjhhGelrStqsVERk0Egn6ntLSL6PN23L3de5e6u6leXl57+Spb9HU1qE59CIigUSCvgoojHtcAFRfRpsB03SuQ5+KFREJJBL0W4ESM5tuZhnAKmBDtzYbgPuD2TfLgIi71/RxrQlrbuskc5Su6EVEIIGgd/dOYA3wArAf+KG77zOzh8zsoaDZRqACKAe+Bfzpheeb2feBV4FZZlZlZh/v4z68RWzoRlf0IiIACV32uvtGYmEef2xt3H0HHr7Ec+95NwVejqa2DmZOGjfQLysiMiiF8pOxTec6tfyBiEggdEHv7jRr6EZE5KLQBf3Z9i6ijt6MFREJhC7oLyxRrOmVIiIxoQv646dbAcjPGpXkSkREBofQBf3h2mYAzboREQmEL+hPtTBuxDAmZY5IdikiIoNC+IK+tpkZk8ZipgXNREQghEFfXtvCzIkathERuSBUQd9wtp36lnZKJo1NdikiIoNGqIL+8KnYG7EzJiroRUQuCFfQ17YAUKIZNyIiF4Uq6MtrWxiTkc6U8SOTXYqIyKARqqA/dKqZGZPGacaNiEicUAX94doWSnunlIkAAASzSURBVDQ+LyLyJgkFvZktN7ODZlZuZo/0cN7M7BvB+d1mtijR5/aVjq4oN5XkccOM3P56CRGRIanXJR7NLB14DLiD2N6wW81sg7u/HtdsBVAS3JYCjwNLE3xunxiensY/fGhBX39bEZEhL5Er+iVAubtXuHs78AywslublcDTHrMZyDKzyQk+V0RE+lEiQZ8PVMY9rgqOJdImkeeKiEg/SiToe5rC4gm2SeS5sW9gttrMysysrK6uLoGyREQkEYkEfRVQGPe4AKhOsE0izwXA3de5e6m7l+bl5SVQloiIJCKRoN8KlJjZdDPLAFYBG7q12QDcH8y+WQZE3L0mweeKiEg/6nXWjbt3mtka4AUgHVjv7vvM7KHg/FpgI3AXUA60Ag++3XP7pSciItIjc+9xyDypSktLvaysLNlliIgMGWa2zd1LezoXqk/GiojIWw3KK3ozqwOOXebTc4H6PixnKFCfU4P6nBout89T3b3HmSyDMujfDTMru9SvL2GlPqcG9Tk19EefNXQjIhJyCnoRkZALY9CvS3YBSaA+pwb1OTX0eZ9DN0YvIiJvFsYrehERiaOgFxEJudAE/UDtZDXQzGy9mdWa2d64Yzlm9jMzOxx8zY4798XgZ3DQzO5MTtXvjpkVmtlLZrbfzPaZ2WeC46Htt5mNNLPXzGxX0OevBsdD2+cLzCzdzHaY2U+Cx6Hus5kdNbM9ZrbTzMqCY/3bZ3cf8jdi6+gcAYqBDGAXMCfZdfVR324CFgF74479HfBIcP8R4G+D+3OCvo8Apgc/k/Rk9+Ey+jwZWBTcHwccCvoW2n4TW9J7bHB/OLAFWBbmPsf1/b8B/wL8JHgc6j4DR4Hcbsf6tc9huaIP7U5W7r4JaOh2eCXwneD+d4Dfjzv+jLufd/c3iC0yt2RACu1D7l7j7tuD+83AfmIb1oS23x7TEjwcHtycEPcZwMwKgLuBJ+IOh7rPl9CvfQ5L0KfaTlaTPLYMNMHXicHx0P0czGwacA2xK9xQ9zsYwtgJ1AI/c/fQ9xn4OvAFIBp3LOx9duBFM9tmZquDY/3a516XKR4iEt7JKuRC9XMws7HAj4DPunuTWU/dizXt4diQ67e7dwELzSwLeM7M5r5N8yHfZzN7P1Dr7tvM7OZEntLDsSHV58D17l5tZhOBn5nZgbdp2yd9DssVfcI7WYXEqWDzdYKvtcHx0PwczGw4sZD/nrv/ODgc+n4DuHsj8DKwnHD3+XrgA2Z2lNhw661m9l3C3WfcvTr4Wgs8R2wopl/7HJagT7WdrDYAHw3ufxT497jjq8xshJlNB0qA15JQ37tisUv3J4H97v6PcadC228zywuu5DGzUcDtwAFC3Gd3/6K7F7j7NGL/Zn/p7vcR4j6b2RgzG3fhPvA+YC/93edkvwPdh+9k30VsdsYR4MvJrqcP+/V9oAboIPa/+8eBCcAvgMPB15y49l8OfgYHgRXJrv8y+3wDsV9PdwM7g9tdYe43MB/YEfR5L/AXwfHQ9rlb/2/md7NuQttnYjMDdwW3fReyqr/7rCUQRERCLixDNyIicgkKehGRkFPQi4iEnIJeRCTkFPQiIiGnoBcRCTkFvYhIyP1/59t/5aMJIDAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# soft constraints\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import cvxpy as cp\n",
    "\n",
    "N = 1\n",
    "T = 500\n",
    "m = 20  # number of devices\n",
    "n = 5  # number of servers\n",
    "y_max = 110\n",
    "y_min = 90\n",
    "rs = 0.1\n",
    "beta = np.ones(n)\n",
    "beta[0] = 3\n",
    "beta[1] = 3\n",
    "gamma = 0.0001\n",
    "\n",
    "def optimization(m, n, wS, y, mu, BS, beta):\n",
    "    x = cp.Variable((n + 1, m), nonneg=True)\n",
    "    obj = beta @ cp.inv_pos(cp.sqrt(x[1:, :] @ y + wS)) + gamma * cp.sum(cp.multiply(cp.diag(mu @ x[1:, :]), y))\n",
    "\n",
    "    constraints = [0 <= x, x <= 1,\n",
    "                   x[1:, :] @ y <= BS,\n",
    "                   cp.sum(x, 0) == 1]\n",
    "\n",
    "    prob = cp.Problem(cp.Minimize(obj), constraints)\n",
    "    prob.solve()  # Returns the optimal value.\n",
    "    return x, prob\n",
    "\n",
    "def oracle(y, mu):\n",
    "    x, prob = optimization(m, n, wS, y, mu, BS, beta)\n",
    "    return x.value, prob.value, prob.status\n",
    "\n",
    "def f(x, y, mu):\n",
    "    return beta.dot(1/np.sqrt(x[1:, :].dot(y) + wS)) + gamma * np.sum((y*x[1:]).T*mu)\n",
    "\n",
    "def f_drop(x, y, mu, BS, gamma):\n",
    "    dp = np.minimum(x[1:, :].dot(y), BS)\n",
    "    return beta.dot(1/np.sqrt(dp + wS)) + gamma * np.sum((y*x[1:]).T*mu)\n",
    "\n",
    "reg = np.zeros((N, T))\n",
    "\n",
    "#records all y, x_opt, x_t #yuhang yao\n",
    "y_N_T = np.zeros((N, T, m)) #yuhang yao\n",
    "x_opt_N_T = np.zeros((N, T, n + 1, m)) #yuhang yao\n",
    "x_t_N_T = np.zeros((N, T, n + 1, m)) #yuhang yao\n",
    "j_N_T = np.zeros((N, T, m))#yuhang yao\n",
    "BS_N = np.zeros((N, n))#yuhang yao\n",
    "\n",
    "for u in range(N):\n",
    "    wS = np.random.randint(15, 25, n)\n",
    "    BS = np.random.uniform(y_min*10, y_max*10, n)\n",
    "    BS_N[u] = BS\n",
    "    mu = np.random.rand(m, n)\n",
    "    mu[:,0] = 0.5\n",
    "    mu[:,1] = 0.5\n",
    "    mu[:,2] = 0.8\n",
    "    mu[:,3] = 0.8\n",
    "    mu[:,4] = 0.8\n",
    "    # trace_gen = Trace(m, n, seed + u)\n",
    "    # mu = trace_gen.avg()\n",
    "    # mu = np.random.rand(m, n)\n",
    "    # mu_hat = np.zeros_like(mu)  # empirical mean\n",
    "    mu_hat = np.zeros_like(mu)\n",
    "    T_ij = np.ones_like(mu)  # total number of times arm (i,j) is played\n",
    "    for t in range(T):\n",
    "        y = np.random.uniform(y_min, y_max, m).astype(int)\n",
    "        x_opt, f_opt, status = oracle(y, mu)\n",
    "        if 'optimal' not in status:\n",
    "            print('Solution infeasible 1')\n",
    "            break\n",
    "\n",
    "        rho_ij = np.sqrt(3 * np.log(t + 1) / (2 * T_ij)) * rs\n",
    "        mu_bar = np.maximum(mu_hat - rho_ij, 0) # LCB\n",
    "        x_t, _, status = oracle(y, mu_bar)\n",
    "        if 'optimal' not in status:\n",
    "            print('Solution infeasible 2')\n",
    "            break\n",
    "\n",
    "        f_t = f(x_t, y, mu)\n",
    "\n",
    "        # sample j based on x_t[i], observe c_ij, update mu_hat[i,j]\n",
    "        # c = trace_gen.generate()\n",
    "        for i in range(m):\n",
    "            j = np.random.choice(n+1, p=x_t[:, i])\n",
    "            j_N_T[u, t, i] = j #yuhang yao\n",
    "            if j != 0:\n",
    "                j -= 1\n",
    "                c_ij = int(np.random.rand() < mu[i, j])\n",
    "                # a = np.random.rand() * 3\n",
    "                # c_ij = np.random.beta(a, a * (1-mu[i, j])/mu[i, j]) # beta distribution\n",
    "                # c_ij = c[i, j]  # trace\n",
    "                T_ij[i, j] += 1\n",
    "                mu_hat[i, j] += (c_ij - mu_hat[i, j]) / T_ij[i, j]\n",
    "\n",
    "        # calculate regert\n",
    "        reg[u, t] = f_t - f_opt\n",
    "        y_N_T[u, t] = y#yuhang yao\n",
    "        x_opt_N_T[u, t] = x_opt#yuhang yao\n",
    "        x_t_N_T[u, t] = x_t#yuhang yao\n",
    "        \n",
    "plt.plot(np.cumsum(reg, axis=1).T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4212202218464382"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta.dot(1/np.sqrt(x_t[1:, :].dot(y) + wS))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([647.11743944, 636.9302228 , 262.27559505, 263.05699937,\n",
       "       268.61974145])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_t[1:, :].dot(y) + wS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([647.11743944, 636.9302228 , 262.27559505, 263.05699937,\n",
       "       268.61974145])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_t[1:, :].dot(y) + wS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1207.3856998158342"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s = 0\n",
    "for i in range(m):\n",
    "    for j in range(n):\n",
    "        s += mu[i,j]*x_t[j+1,i]*y[i]\n",
    "s"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
